{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Self-organized criticality\n",
    "\n",
    "Code examples from [Think Complexity, 2nd edition](http://greenteapress.com/wp/complexity2), Chapter 8\n",
    "\n",
    "Copyright 2016 Allen Downey, [MIT License](http://opensource.org/licenses/MIT)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "from __future__ import print_function, division\n",
    "\n",
    "%matplotlib inline\n",
    "\n",
    "import numpy as np\n",
    "import matplotlib.pyplot as plt\n",
    "\n",
    "import thinkplot\n",
    "from thinkstats2 import Hist, Cdf\n",
    "from thinkstats2 import RandomSeed\n",
    "\n",
    "from matplotlib import rc\n",
    "rc('animation', html='html5')"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Sand pile"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "`Sand.py` contains an implementation of the sand pile model."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "from Sand import SandPile, SandPileViewer"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Here's a small example starting with two cells ready to topple.  `n` is the number of rows, `m` is the number of columns."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "pile = SandPile(n=3, m=5, level=0)\n",
    "pile.array[1, 1] = 4\n",
    "pile.array[1, 3] = 4\n",
    "\n",
    "a = pile.array\n",
    "print(a)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "To execute one step, first we find cells that are above the toppling threshold, `K`."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "K = 3\n",
    "toppling = a > K\n",
    "print(toppling.astype(int))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Then we use `correlate2d` to make a copy of the update kernel around each toppling cell."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "kernel = np.array([[0, 1, 0],\n",
    "                   [1,-4, 1],\n",
    "                   [0, 1, 0]], \n",
    "                  dtype=np.int32)\n",
    "\n",
    "from scipy.signal import correlate2d\n",
    "\n",
    "c = correlate2d(toppling, kernel, mode='same', boundary='fill', fillvalue=0)\n",
    "print(c)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Finally, we add the result back into the array:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "a += c\n",
    "print(a)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Any grains that topple off the edge disappear."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Animation"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Let's look at a bigger pile, with `n=20`.  All cells are initialized to `level`, which is meant to be substantially bigger than `K`."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "pile = SandPile(n=20, level=10)\n",
    "print(pile.run())"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "The `run` function invokes `step` until no more cells topple and returns the number of time steps and the number of affected cells."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Now let's look at an animation, starting from this initialized pile."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "viewer = SandPileViewer(pile, drop_flag=True)\n",
    "anim = viewer.animate(frames=100)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Starting with `level>K` produces all kinds of interesting patterns."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Each step of the animation drops a single grain at a random location and runs until no more cells topple."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Most avalanches are small, but some are very large."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "anim"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "After a while, the pile looks pretty random.\n",
    "\n",
    "Here's a plot of the number of cells toppled after each `step`.  "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "thinkplot.plot(pile.toppled_seq)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "The following figure shows the progression of the pile from ordered to apparently random."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "RandomSeed(17)\n",
    "\n",
    "pile = SandPile(n=20, level=10)\n",
    "print(pile.run())\n",
    "viewer = SandPileViewer(pile)\n",
    "\n",
    "thinkplot.preplot(cols=3)\n",
    "viewer.draw()\n",
    "\n",
    "thinkplot.subplot(2)\n",
    "for i in range(200):\n",
    "    viewer.step()\n",
    "viewer.draw()\n",
    "\n",
    "thinkplot.subplot(3)\n",
    "for i in range(200):\n",
    "    viewer.step()\n",
    "viewer.draw()\n",
    "\n",
    "plt.savefig('chap08-1.pdf')"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Long tailed distributions"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "If the sand pile is in a critical state, we expect quantities like the duration of an avalanche, `T`, and the number of cells affected, `S`, to have long-tailed distributions.\n",
    "\n",
    "Following Bak, Tang, and Wiesenseld, we start with a 50x50 array and plot the PMFs of `S` and `T` on a log-log scale."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "scrolled": true
   },
   "outputs": [],
   "source": [
    "pile2 = SandPile(n=50, level=30)\n",
    "pile2.run()\n",
    "viewer = SandPileViewer(pile2)\n",
    "viewer.draw()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Now we run the pile for many time steps and keep track of the durations and number of cells affected."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "scrolled": true
   },
   "outputs": [],
   "source": [
    "RandomSeed(17)\n",
    "\n",
    "iters = 100000\n",
    "%time res = [pile2.drop_and_run() for _ in range(iters)]"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "And extract the results into two NumPy arrays."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "T, S = np.transpose(res)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "A large majority of drops have duration 1 and no toppled cells, so we'll filter those out."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "T = T[T>1]\n",
    "S = S[S>0]"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "The distributions of `T` and `S` have many small values and a few very large ones.\n",
    "\n",
    "Here are the histograms on linear axes, showing only small values."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "thinkplot.preplot(cols=2)\n",
    "\n",
    "histT = Hist(T)\n",
    "thinkplot.Pdf(histT, label='T')\n",
    "thinkplot.config(xlabel='Avalanche duration',\n",
    "                 ylabel='Number',\n",
    "                 xlim=[1, 50], loc='upper right')\n",
    "\n",
    "thinkplot.subplot(2)\n",
    "\n",
    "histS = Hist(S)\n",
    "thinkplot.Pdf(histS, label='S')\n",
    "thinkplot.config(xlabel='Avalanche size',\n",
    "                 ylabel='Number',\n",
    "                 xlim=[1, 50])\n",
    "\n",
    "plt.savefig('chap08-2.pdf')"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "To see whether these distributions follow a power law, we plot the histograms on a log-log scale."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "thinkplot.preplot(cols=2)\n",
    "\n",
    "x = 9000\n",
    "thinkplot.plot([1, x], [x, 1], color='gray', linewidth=3)\n",
    "thinkplot.Pdf(histT, label='T', linewidth=2)\n",
    "thinkplot.config(xlabel='Avalanche duration',\n",
    "                 xlim=[1, 1000],\n",
    "                 ylabel='Number',\n",
    "                 xscale='log',\n",
    "                 yscale='log',\n",
    "                 loc='upper right')\n",
    "\n",
    "thinkplot.subplot(2)\n",
    "\n",
    "x = 5700\n",
    "thinkplot.plot([1, x], [x, 1], color='gray', linewidth=3)\n",
    "thinkplot.Pdf(histS, label='S', linewidth=1)\n",
    "thinkplot.config(xlabel='Avalanche size',\n",
    "                 xlim=[1, 5600],\n",
    "                 ylabel='Number',\n",
    "                 xscale='log',\n",
    "                 yscale='log')\n",
    "\n",
    "plt.savefig('chap08-3.pdf')"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "The gray line has slope -1.  The distribution of avalanche duration is approximately straight between 1 and 100, but then drops off.  The distribution of size follows a power law more closely and over a greater range, but it also seems to drop off for values above a few hundred."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "**Exercise:** Try running the model longer to see if you can get a less noisy plot of the distributions of `T` and `S`."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Fractals\n",
    "\n",
    "If the sand pile is in a critical state, we expect to see fractal geometry.\n",
    "\n",
    "To estimate the fractal dimension, I'll start with a bigger pile and a higher initial level."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "pile3 = SandPile(n=131, level=22)\n",
    "%time pile3.run()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "The initial state sure looks like a fractal."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "viewer3 = SandPileViewer(pile3)\n",
    "viewer3.draw()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Since it contains four different levels (0, 1, 2, and 3), we can extract 4 binary patterns."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "def draw_four(viewer, vals=range(4)):\n",
    "    thinkplot.preplot(rows=2, cols=2)\n",
    "    a = viewer.viewee.array\n",
    "    \n",
    "    for i, val in enumerate(vals):\n",
    "        thinkplot.subplot(i+1)\n",
    "        viewer.draw_array(a==vals[i], vmax=1)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Here's what they look like:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "draw_four(viewer3)\n",
    "plt.savefig('chap08-4.pdf')"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Now we can apply a box-counting algorithm to each level.\n",
    "\n",
    "`count_cells` starts with a single cell in the middle, gradually increases the size of the box, and counts the number of cells in each box."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "def count_cells(a):\n",
    "    \"\"\"Counts the number of cells in boxes with increasing size.\n",
    "    \n",
    "    a: NumPy array\n",
    "    \n",
    "    returns: list of (i, i**2, cell count) tuples\n",
    "    \"\"\"\n",
    "    n, m = a.shape\n",
    "    end = min(n, m)\n",
    "    \n",
    "    res = []\n",
    "    for i in range(1, end, 2):\n",
    "        top = (n-i) // 2\n",
    "        left = (m-i) // 2\n",
    "        box = a[top:top+i, left:left+i]\n",
    "        total = np.sum(box)\n",
    "        res.append((i, i**2, total))\n",
    "        \n",
    "    return np.transpose(res)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "`box_count` takes a pile and a value, extracts the cells that have the given value, calls `count_cells`, and estimates the fractal dimension.\n",
    "\n",
    "If `plot` is `True`, it also generates a graph of cell count versus box size on a log-log scale."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "from scipy.stats import linregress\n",
    "\n",
    "def box_count(pile, val=0, plot=False):\n",
    "    \"\"\"Estimates the fractal dimension by box counting.\n",
    "    \n",
    "    pile: SandPile\n",
    "    val: which value from the pile to count\n",
    "    plot: boolean, whether to generate a plot\n",
    "    \n",
    "    returns: estimated fractal dimension\n",
    "    \"\"\"\n",
    "    res = count_cells(pile.array==val)\n",
    "    steps, steps2, cells = res\n",
    "    \n",
    "    # select the range where we have a nonzero number of cells\n",
    "    legit = np.nonzero(cells)\n",
    "    steps = steps[legit]\n",
    "    steps2 = steps2[legit]\n",
    "    cells = cells[legit]\n",
    "\n",
    "    if plot:\n",
    "        thinkplot.preplot(3)\n",
    "        thinkplot.plot(steps, steps2, linestyle='dashed')\n",
    "        thinkplot.plot(steps, cells, label='val=%d' % val)\n",
    "        thinkplot.plot(steps, steps, linestyle='dashed')\n",
    "\n",
    "        thinkplot.config(xscale='log', yscale='log',\n",
    "                         xlim=[1, 200], loc='upper left',\n",
    "                         xlabel='box size', ylabel='cell count')\n",
    "    \n",
    "    params = linregress(np.log(steps), np.log(cells))\n",
    "    return params[0]"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Finally `box_count_four` applies the box counting algorithm for each value in the sand pile."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "def box_count_four(pile, vals=range(4)):\n",
    "    \"\"\"Applies box counting to each value in the pile.\n",
    "    \n",
    "    pile: SandPile\n",
    "    vals: list of values to check\n",
    "    \"\"\"\n",
    "    thinkplot.preplot(rows=2, cols=2)\n",
    "    \n",
    "    dims = []\n",
    "    for i, val in enumerate(vals):\n",
    "        thinkplot.subplot(i+1)\n",
    "        dim = box_count(pile, val, plot=True)\n",
    "        dims.append(dim)\n",
    "        \n",
    "    return dims"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Here are the results:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "dims = box_count_four(pile3)\n",
    "plt.savefig('chap08-5.pdf')"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "The lines are reasonably straight, which indicates that we are running the algorithm over a valid range of box sizes.  Here are the estimated slopes:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "for i, dim in enumerate(dims):\n",
    "    print('%d  %0.3f' % (i, dim))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "The fractal dimension for values 0, 1, and 2 seems to be clearly non-integer, which indicates that the image is fractal.\n",
    "\n",
    "Strictly, the fractal dimension for value 3 is indistinguishable from 2, but given the results for the other values, the apparent curvature of the line, and the appearance of the pattern, it seems likely that it is also fractal. "
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "**Exercise:**  Choose a different value of `n` and/or `level` and run this analysis again.  Are the estimated fractal dimensions consistent?"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Spectral density"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Suppose the sandpile made a little click each time a cell toppled.  What would it sound like?\n",
    "\n",
    "`toppled_seq` contains the number of cells that toppled during each time step.  We can use Welch's algorithm to estimate its [power spectral density](https://en.wikipedia.org/wiki/Spectral_density)."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "signal = pile2.toppled_seq\n",
    "len(signal)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "from scipy.signal import welch\n",
    "\n",
    "nperseg = 2048\n",
    "freqs, powers = welch(signal, nperseg=nperseg, fs=nperseg)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "And here's what it looks like."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "x = nperseg\n",
    "ys = np.array([x**1.58, 1]) / 2.7e3\n",
    "thinkplot.plot([1, x], ys, color='gray', linewidth=1)\n",
    "\n",
    "thinkplot.plot(freqs, powers)\n",
    "thinkplot.config(xlabel='frequency', xscale='log', xlim=[1, 1200],\n",
    "                 ylabel='power', yscale='log', ylim=[1e-5, 5])\n",
    "\n",
    "plt.savefig('chap08-6.pdf')"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "The slope of the line is -1.58, which indicates that this spectrum is pink noise with parameter $\\beta=1.58$."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "**Exercise:**  Choose a different value of `nperseg` and run this analysis again.  What are the pros and cons of larger segment lengths?  Modify the code to run the model longer and see if you can get a less noisy estimate of the spectrum."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Exercises"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "**Exercise:** To test whether the distributions of `T` and `S` are heavy-tailed, we plotted their histograms on a log-log scale, which is what Bak, Tang and Wiesenfeld show in their paper.  But as we saw in Chapter 4, this visualization can obscure the shape of the distribution.  Using the same data, make a plot that shows the CDFs of `S` and `T`.  What can you say about the shape of these distributions?  Do they follow a power law?  Are they heavy tailed?\n",
    "\n",
    "You might find it helpful to plot the CDFs on a log-x scale and on a log-log scale."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Solution\n",
    "\n",
    "# Here is the CDF of `S` on a log-x scale (left) and on a log-log scale (right).\n",
    "\n",
    "thinkplot.preplot(cols=2)\n",
    "\n",
    "cdf = Cdf(S)\n",
    "thinkplot.Cdf(cdf)\n",
    "thinkplot.config(xlabel='Avalanche size',\n",
    "                 ylabel='Cdf',\n",
    "                 xscale='log')\n",
    "\n",
    "thinkplot.subplot(2)\n",
    "                  \n",
    "thinkplot.Cdf(cdf, complement=True)\n",
    "thinkplot.config(xlabel='Avalanche size',\n",
    "                 ylabel='Cdf',\n",
    "                 xscale='log',\n",
    "                 yscale='log')"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Solution\n",
    "\n",
    "# And here are the results for `T`.\n",
    "\n",
    "thinkplot.preplot(cols=2)\n",
    "\n",
    "cdf = Cdf(T)\n",
    "thinkplot.Cdf(cdf)\n",
    "thinkplot.config(xlabel='Avalanche size',\n",
    "                 ylabel='Cdf',\n",
    "                 xscale='log')\n",
    "\n",
    "thinkplot.subplot(2)\n",
    "                  \n",
    "thinkplot.Cdf(cdf, complement=True)\n",
    "thinkplot.config(xlabel='Avalanche size',\n",
    "                 ylabel='Cdf',\n",
    "                 xscale='log',\n",
    "                 yscale='log')"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Solution\n",
    "\n",
    "# On a log-log scale, these CDFs are not straight, \n",
    "# which indicates that they do not follow a power law.\n",
    "\n",
    "# However, on a log-x scale, they resemble lognormal distributions\n",
    "# (with the left tail truncated).  The lognormal distribution is\n",
    "# also considered heavy-tailed.\n",
    "\n",
    "# In my opinion, many papers in complexity science are too quick to \n",
    "# identify power laws based on plots showing histograms or probability \n",
    "# mass functions (PMFs).  I think many of these distributions, \n",
    "# including both measurements from the world and results from simulations,\n",
    "# are better modeled by other heavy-tailed distributions.\n",
    "\n",
    "# However, the argument the sand pile model is intended to make is valid \n",
    "# either way: we observe that heavy-tailed distributions are common in \n",
    "# natural systems, and SOC is a possible explanation, regardless of whether\n",
    "# the distributions are specifically power law distributions or more generally \n",
    "# heavy-tailed."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "**Exercise:** In Section 8.5 we showed that the initial configuration of the sand pile model produces fractal patterns.  But after we drop a large number of random grains, the patterns look more random.\n",
    "\n",
    "Starting with the example in Section 8.5, run the sand pile model for a while and then compute fractal dimensions for each of the 4 levels.  Is the sand pile model fractal in steady state?"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Solution\n",
    "\n",
    "# I'll start with a copy of `pile4`.  Again the starting configuration looks fractal.\n",
    "\n",
    "from copy import deepcopy\n",
    "\n",
    "pile4 = deepcopy(pile3)\n",
    "viewer4 = SandPileViewer(pile4)\n",
    "draw_four(viewer4)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Solution\n",
    "\n",
    "# After it runs for a while, it looks pretty random.\n",
    "# But we have to be careful; people are not very good at judging \n",
    "# what is random or not.\n",
    "\n",
    "for i in range(10000):\n",
    "    pile4.drop_and_run()\n",
    "    \n",
    "draw_four(viewer4)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Solution\n",
    "\n",
    "# Here are the results from the box counting algorithm.\n",
    "# The lines are reasonably straight, which means our\n",
    "# estimates of the fractal dimension are valid.\n",
    "\n",
    "dims = box_count_four(pile4)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Solution\n",
    "\n",
    "# And here are the fractal dimensions.\n",
    "\n",
    "for i, dim in enumerate(dims):\n",
    "    print('%d  %0.3f' % (i, dim))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Solution\n",
    "\n",
    "# They all seem to be converging on 2, which is what we expect\n",
    "# from a totally random pattern.  It's hard to say for sure, \n",
    "# because it's possible that they are still fractal with dimensions \n",
    "# close to 2, but not exactly 2.  But I think that's unlikely.\n",
    "\n",
    "# I think this observation undermines the SOC argument to some degree.\n",
    "# If the sand pile model is not fractal in steady state, it does not \n",
    "# explain the prevalence of fractal geometry in natural systems."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "**Exercise:** Another version of the sand pile model, called the \"single source\" model, starts from a different initial condition: instead of all cells at the same level, all cells are set to 0 except the center cell, which is set to a very large value.\n",
    "\n",
    "Write a function that creates a `SandPile` object, sets up the single source initial condition, and runs until the pile reaches equilibrium.  Does the results appear to be fractal?\n",
    "\n",
    "You can read more about this version of the sand pile model at http://math.cmu.edu/~wes/sandgallery.html"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Solution\n",
    "\n",
    "# The following function clears all cells and puts a tower in the center.\n",
    "    \n",
    "def single_source(pile, height=1024):\n",
    "    \"\"\"Adds a tower to the center cell.\n",
    "    \n",
    "    height: value assigned to the center cell\n",
    "    \"\"\"\n",
    "    a = pile.array\n",
    "    n, m = a.shape\n",
    "    a[:, :] = 0\n",
    "    a[n//2, m//2] = height"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Solution\n",
    "\n",
    "# `run_ss_pile` builds a single-source sand pile and runs until equilibrium.\n",
    "\n",
    "def run_ss_pile(n, k=10):\n",
    "    \"\"\"Runs a single source model.\n",
    "    \n",
    "    k: integer power of two of height\n",
    "    \"\"\"\n",
    "    pile = SandPile(n)\n",
    "    single_source(pile, 2**k)\n",
    "    print(pile.run())        \n",
    "    viewer = SandPileViewer(pile)\n",
    "    return pile, viewer"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Solution\n",
    "\n",
    "# Let's try it with `2**14` grains:\n",
    "\n",
    "% time pile5, viewer5 = run_ss_pile(n=101, k=14)\n",
    "viewer5.draw()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Solution\n",
    "\n",
    "# Visually, the results sure look fractal.\n",
    "\n",
    "draw_four(viewer5)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Solution\n",
    "\n",
    "# And here are the results from the box counting algorithm.\n",
    "\n",
    "# For value 0, the line gets steeper at the end because of the boundary cells.\n",
    "\n",
    "# For the other values, it drops off at the end.\n",
    "\n",
    "# So the estimated dimensions might not be very accurate.\n",
    "# A better alternative would be a radial version of the box counting algorithm.\n",
    "\n",
    "dims = box_count_four(pile5)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Solution\n",
    "\n",
    "# The results are not as clear as we might like, \n",
    "# but at least some of these patterns seem to have non-integer dimensions.\n",
    "\n",
    "for i, dim in enumerate(dims):\n",
    "    print('%d  %0.3f' % (i, dim))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Solution\n",
    "\n",
    "# In the repository for this book, you'll find a notebook called `sspile` \n",
    "# that does a more careful job of estimating these dimensions."
   ]
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "Python 3",
   "language": "python",
   "name": "python3"
  },
  "language_info": {
   "codemirror_mode": {
    "name": "ipython",
    "version": 3
   },
   "file_extension": ".py",
   "mimetype": "text/x-python",
   "name": "python",
   "nbconvert_exporter": "python",
   "pygments_lexer": "ipython3",
   "version": "3.6.4"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 1
}
